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Solving coupled nonlinear sine-Gordon equations and Maxwell equations numerically, we study 
the electromagnetic and superconducting properties of the single crystal of high-T c superconductor 
Bi2Sr2CaCu2 08+5 with a static magnetic field applied parallel to the ab- plane and a dc current fed 
in along the c-axis. Cavity resonances of transverse plasma occur in the intrinsic Josephson junctions 
with frequencies in terahertz regime. It is revealed that the electromagnetic wave can transmit from 
the junctions into space. The emitted energy counted by the Poynting vector is about 400W/cm 2 . 
The frequency as well as the energy of emission can be tuned almost continuously by the current 
and magnetic field. 

PACS numbers: 74.50. 74.25.Gz 85.25.Cp 



Introduction - Terahertz (THz) technology is an ex- 
tremely attractive field. The main users of the THz elec- 
tromagnetic waves are perhaps the biomedical diagnos- 
tics, DNA probe and cancer detection; a THz tomogra- 
phy is also very useful for material characterization 
Although the lower and higher frequency bands of elec- 
tromagnetic field can be generated by electronics and 
photonics respectively, seeking solid-state, stable genera- 
tors for the THz waves is still a subject of scientific effort. 
While recently the semiconductor heterostructures gen- 
erate THz emissions with high efficiency Q, known as the 
quantum cascade lasers, the issue of frequency tunability 
has not been resolved. 

In the present work, following previous proposal^, 
we seek continuously tunable emissions of THz elec- 
tromagnetic waves from high-T c superconductors with 
the principle of Josephson relation [J]. For this pur- 
pose, the highly anisotropic high-T c superconductor 
B^S^CaC^Os+a can be considered as a stack of 
Josephson junctions in atomic scale[5|. The novel de- 
vice has its advantages since, first, the energy gap in 
high-r c superconductors is much larger than the plasma 
energy and thus the plasma, if excited in some way, 
should be stable; secondly, the power output conjectured 
to be proportional to the junction number squared would 
be very large; thirdly, variations of physical parameters 
in these intrinsic Josephson junctions (IJJs) are much 
smaller than artificially fabricated junctions. For artifi- 
cial Josephson junction arrays, radiations of electromag- 
netric waves have been demonstrated Q; the frequency 
is, however, in the sub-THz regime because of the small 
superconducting energy gap. 

There are theoretical calculations as well as numeri- 
cal simulations [1, 0] which discussed possible radiation 
from IJJs of high-Xc cuprates. However, many issues 
have not yet been revealed concerning the mechanism 
of THz emission. Although the Josephson plasma ob- 



viously play a key role here, the physical parameters of 
single Josephson junctions of I^S^CaCu^Os+a only give 
plasma frequency in sub-THz regime. Therefore the col- 
lective modes in the stack of Josephson junctions are es- 
sential in order to lift the frequency by an order of magni- 
tude. This goal is expected to be achieved by the Joseph- 
son vortices driven by the c-axis current, which generates 
oscillating voltage across junctions according to the ac 
Josephson relation. The Josephson plasma and the mo- 
tion of vortices intervene with each other in a complex 
way, which makes the resonance condition in this system 
obscure. Without revealing the resonance mechanism, it 
is hard to tell whether a continuous tuning of frequency 
is possible or not. 

On the other hand, experimental efforts toward ex- 
citing THz electromagnetic wave using IJJs seem to be 
accelerated recently [lfj, [lj], [l^]. However, there are 



large discrepancies in estimates of the optimal power out- 
put: theoretically 3000W/cm 2 in Ref. @, 1500W/cm 2 
in Ref. Q, while experimentally 150W/cm 2 in Ref. [Til ] . 
20W/cm 2 in Ref. [12j, which hinders a clear assessment 
on the new technique. 

In this paper, we investigate the THz radiation from 
intrinsic Josephson junctions by solving coupled nonlin- 
ear sine-Gordon equations numerically using an appropri- 
ate boundary condition. The main results of the present 
work are summarized as follows: Resonant radiation of 
THz wave occurs at the edge of junctions due to the cav- 
ity resonance of transverse plasma. A large resonance is 
achieved when the velocity of Josephson vortices matches 
the velocity of Josephson plasma. The vortex configura- 
tion is revealed to be rectangular with additional ran- 
dom sliding motions at resonance. The maximum energy 
is about 400W/cm 2 and the frequency covers the THz 
band. Both the energy and frequency can be tuned al- 
most continuously by the bias current and applied mag- 
netic field. 
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FIG. 1: (color online). Schematic diagram of terahertz elec- 
tromagnetic wave emitter. The intrinsic Josephson junctions 
are sandwiched by two gold electrodes through which a dc 
bias current is fed into the system. 

Model equations - The model we use, sketched schemat- 
ically in Fig. [T] is the same as Ref. 0] . A static magnetic 
field B a of order of IT is applied along the y axis, which 
induces Josephson vortices in all insulating layers. A bias 
current is fed into the system along the z axis, which 
drives fluxons towards the negative direction of the x 
axis. The moving fluxons excite the transverse Joseph- 
son plasma [13l. Il4j . For simplicity, we ignore here the 
effect of thermal fluctuations and thus the vortex lines 
are straight. This approximation makes one to be able 
to reduce the three dimensional system to two dimen- 
sions. While the estimate on the output energy should 
be considered as the upper limit for experiments, the 
mechanism of radiation remains unchanged even when 
thermal fluctuations are involved. 

Using the London theory, Josephson relation as well 
as Maxwell equations, the superconducting and electro- 
magnetic properties of the system can be described by 
the following equations (il flil : 

(1 -(A^){d 2 ,P l+u + pdfPi +u +sinP i+1)i m 
+as'[d t , (p' l+1 - p\) + m+i Pi)]} = d%Pi+i,U { ) 

s'(l - aA^)p[ = d t ,(P l+u - P w _i), (2) 

where Pi+i^i is the guage invariant phase difference de- 
fined as 

Pl +1 l(r,t) = ipi +1 (r,t) - (pi(r,t) - — / A z (r,t)dz, 

<Po J Zl 

(3) 

with the phase of order parameter ip, the vector poten- 
tial A z and the quantum flux <pQ. The operator A^ 2 ) is 
defined as A^ft = f l+1 +/,_! - 2f u In Eqs.Q and ©, 
dimensionless quantities are used: x' — x/X c , t' = ttu p 
and p[ = pi\ c Lo p / J c , where pi is the charge density in the 
Zth superconducting layer, lu p — c/ \ c ^fe2 is the plasma 
frequency, e c is the dielectric constant along the z axis. 



A c , and A a fc in below, are the penetration depths. The 
other dimensionless quantities are defined as: J' = J / J c , 
B' = 2ttA c ZAB/0 o , (3 = 4ttct c A c / ' cJT c , E' = a c E/f3J c , 
a = e c p 2 /sD (capacitive coupling), £ = X 2 ab /sD (induc- 
tive coupling) and s' — s/A c , where a c is the conduc- 
tivity, p is the Debye screening length, J c is the critical 
current density, s (D) is the thickness of the supercon- 
ducting (insulating) layer. 

Boundary condition - Equations |T]) and ((2|) must be 
modified at the topmost and bottommost junctions. We 
assume that the superconductivity will penetrate into 
the gold electrodes [9( so that the thickness of the top- 
most and bottommost superconducting layers so is larger 
than s. We also assume the electrodes as good con- 
ductor, therefore the electric field in the electrodes is 
zero. With these two assumptions, the equations at 
the topmost and bottommost layer can be written down 
straightforwardly (l7j . 

The boundary condition at the edge (x = and 
x = L) is more subtle. We implemented the dynamic 
boundary condition which is determined by the electro- 
magnetic wave in the dielectric medium coupled to the 
junctions^, EH- Assuming that only transverse waves 
are transmitted through the dielectric medium, which 
becomes exact when the number of junctions is infinite, 
there is a linear relation between the electric and mag- 
netic fields. With the Fresnel conditions, one arrives at 
following condition between the scillating fields B' and 
E' on the IJJs side of interface [3| 

B' v = T^ d E ,z , (4) 

where + (— ) means the wave propagating in the nega- 
tive(positive) x direction, e' d = ed/s c is the dielectric con- 
stant of the dielectric medium ea normalized by e c . With 
the relations d x >Pi +1 j = B' v l+l l and dfPi+i.i = E' Z 1+1U 
we derive from Eq.(j4]) the dynamic boundary condition 
for the IJJs 

d x ,P l+U - B' a ± yp^&vPl+Xt - E' dc ) = 0, (5) 

where the magnetic field induced by J^ xt is ignored be- 
cause it is negligibly small in comparison to B' a . In high- 
bias region E' dc w J' cxt j 1 (3 (see [9(). 

The parameters used for simulations are \ a b = 0.4^m, 
A c = 200Aim, s = 3A, s = 0.075/jm, D = 12A, p = 0.6A, 
a = 0.1, (3 = 0.02 and s d = e c = 10. The number 
of junctions is N = 30 and the length is L = 19.2^m. 
The applied magnetic field is B a = IT. The equations of 
motion Eqs.([T]) and (J2) are integrated by the five-point 
Nordsieck-Gear Predictor-corrector method. The time 
step in all the simulations is set to At 1 = 0.0002 and 
the mesh size is set to Ax' = 0.0002 (0.1A a b). We start 
to calculate physical quantities when the system reaches 
steady states which are identified by the criterion that 
the voltage only fluctuates weakly around a fixed value. 
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FIG. 2: Radiation power counted by the Poynting vector as a 
function of the voltage across the stack of IJJs. The external 
current corresponds to the b-peak is J ex t = 1.33J C . The inset 
is the enlarged view of the main panel in the region V < 
235mV. 



Numerical results -We gradually ramp up the current 
J' cy± at B a = IT and then gradually reduce it with step 
AJg Xt =0.01. The output power is counted by the Poynt- 
ing vector at the left edge, toward which the vortices are 
moving. The dependence of power on voltage V is plotted 
in Fig. [2] The resonance takes place at discrete values 
of voltage with equidistance between neighboring peaks. 
The maximum power is about 400W/cm 2 . 

The spectra of the electric field at the main peaks in 
Fig. [2] obtained by fast Fourier transformation are dis- 
played in Fig. [3l The sharp spectra indicate monochro- 
matic electromagnetic waves in terahertz regime. Thus 
we have confirmed theoretically the terahertz laser radi- 
ation from the intrinsic Josephson junctions. 

The equidistance between neighboring peaks indicates 
a resonance phenomenon, the one between Josephson 
plasma and cavity modes as revealed below. The Joseph- 
son relation reads wjp = 2eV/h and the cavity modes 
read w c = mrc p /L where c p is the velocity of transverse 
plasma and n is an integer. The resonance occurs when 
wjp = oj Ci namely 



V = V/N = nirhc p /2eL. 



(6) 



To be more specific, we counted the wave length A^, fre- 
quency / and the voltage V for the resonating peaks 
in Fig. [§J a-peak: X w = 1.82/mi, / = 4.79THz, V = 
9.90mV; b-peak: X w = 1.74/mi, / = 5.01THz, V = 
10.36mV and c-peak: X w = 1.66/xm, / = 5_.24THz, 
V = 10.83mV. It is easy to see that / = 2eV/h and 
nX w /2 = L for all the cases: a-peak, n = 21; b-peak, 
n = 22; c-peak, n = 23. The former is nothing but the 
Josephson relation and the latter is the condition for cav- 



ity modes. From the wave length and the frequency, we 
estimate the plasma velocity as c p w 8.71 x 10 6 m/s. It 
is then easy to see that the interval between consecutive 
peaks AV = irhc p /2eL = 0.47mV is satisfied perfectly 
Therefore our simulation indicates clearly that the tera- 
hertz laser is caused by the cavity resonance. 

The above plasma velocities c p is very close to the 
largest characteristic velocity c\ — 8.75 x 10 6 m/s ob- 
tained by solving the linearized Eqs.([T]) and (J2J) with 
so > s (see also [19(). The velocity c\ associated with the 
node-less mode is most important for resonance since it 
corresponds to the in-phase motion of vortices. 

From Fig. [3J it is clear that there is a bundle of 
large resonances around V — 310. 8mV for B a = IT. 
In order to reveal the reason for this nonlinear prop- 
erty, we investigate the vortex motion. Since P/+i,f = 
2wDB a x/(j)o + 2ncDE[ +l l t/(j) + Pi+ij, where Pi+i,i is 
the oscillating contribution, the velocity of vortices is 



cEf +ll /B a = 8.63 x 10 fa m/s at the 



evaluated as v 
b-peak in Fig. [2l It is close to the plasma velocity c p and 
C\. Therefore, the largest resonance takes place when 
the velocities of vortices and transverse plasma coincide 
with each other. A similar relation has been discussed in 
Ref. [2Cj| . where an infinitely long junction was considered 
and therefore no cavity mode was involved. 

From v — c p = C\, it becomes clear that the largest 
energy emission is excited by a voltage 



V = N Cl B a D/c. 



(7) 



With Eqs. © and Q, we have n = 2B a LD/(j> . Since 
nX w /2 = L, it is concluded that the optimal output is 
achieved when the wave length of plasma equals to the 
vortex- vortex separation A„, = <po/B a D. 
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FIG. 3: (color online). Frequency spectra for the voltage at 
the main peaks in Fig. [2] 
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FIG. 4: (color online). Left panel: snapshots for the magnetic 
and electric fields at the largest resonance in Fig. [5] where 
the static backgrounds have been subtracted. Right panel: 
snapshot of vortex configuration. 



It is interesting to take a look at the vortex configura- 
tion at resonance at this point. As seen from the snapshot 
shown in Fig. ([4]), the vortices form an overall ordered 
rectangular lattice; sliding motions take place in random 
places for the time being. The rectangular vortex lattice 
is in accordance with the matching between the measured 
electromagnetic wave velocity and the node-less plasma 
velocity. 

Having clarified the mechanism of resonance, we turn 
to investigate how to tune the resonance frequency and 
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power. The voltage dependence of power emission at dif- 
ferent magnetic fields is presented in Fig. [5] The voltage 
for largest power output increases linearly with the mag- 
netic field, consistent with Eq. [7] From the ac Josephson 
relation, as depicted in the upper axis of Fig. [51 the fre- 
quency of emitted electromagnetic wave can be tuned by 
sweeping the magnetic field and adjusting the dc volt- 
age accordingly. Although the cavity resonance imposes 
that the resonating frequency is discrete when voltage 
is swept, it becomes almost continuous if the length of 
IJJs L is large enough. For example, the frequency in- 
terval between neighboring resonances is about 0.04THz 
for L — 100/im. 

At a given voltage, on the other hand, the emitted 
power shrinks quickly when the magnetic field is tuned 
away from the optimal matching value. This permits one 
to control the power output in a sensitive way. 

Discussions - Finally we compare briefly our results with 



160 240 320 400 480 
Voltage [mV] 



recent experimental works. Both Refs. and [12 1 
tried to stimulate THz electromagnetic waves utilizing 
the Josephson vortex dynamics. While the emission was 
detected by bolometer in Ref. 11], Ref. 12] used a sec- 
ond stack of IJJs to detect excited THz electromagnetic 
waves. It is somewhat difficult to draw conclusions on 
the resonance mechanism from experiments; in our sim- 
ulation, however, cavity resonances are clearly observed. 
In Ref. [jjj it is found that the most efficient emission is 
achieved for the rectangular vortex configuration, same 
as our computer simulation. Although not revealing a 
clear relation, both works tried to tune the frequency 
of emission by sweeping simultaneously the voltage and 
magnetic field, which, according to our theoretical inves- 
tigation, is an efficient way for achieving the frequency 
tunability. In our simulation, we use magnetic fields 
below 1.5 Tesla, which are quite smaller than those in 
Ref. [13]. According to our relation for the optimal 
matching magnetic field and voltage, in order to get large 
emission power under 4 Tesla the voltage should be un- 
rcalistically large. It can be a possible reason that our 
theoretical estimate on the power output is different from 
that in the experiment. 
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FIG. 5: (color online). Radiation power at different applied 
magnetic fields in units of Telsa as a function of the voltage 
across the stack of IJJs. The frequency in the upper axis is 
determined by the ac Josephson relation. 
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